Liquid-liquid transition in supercooled silicon determined by 

first-principles simulation 



P. Ganesh 

Carnegie Institution of Washington, Washington DC 
M. Widom 

Department of Physics, Carnegie Mellon University, Pittsburgh, PA 
(Dated: September 16, 2008) 

First principles molecular dynamics simulations reveal a liquid-liquid phase tran- 
sition in supercooled elemental silicon. Two phases coexist below Tc ~ 1232i^. 
The low density phase is nearly tetra-coordinated, with a pseudogap at the Fermi 
surface, while the high density phase is more highly coordinated and metallic in 
nature. The transition is observed through the formation of van der Waals loops in 
pressure- volume isotherms below Tc. 

PACS numbers: 

Silicon occupies a position in the periodic table at the border between metals and insu- 
lators. At low pressure, crystalline silicon (c-Si) is tetrahedrally coordinated (like diamond) 
and is an indirect band-gap semiconductor. As pressure increases, crystalline Si trans- 
forms from the diamond phase through a more highly-coordinated metallic /5-tin phase to 
a hexagonal phase. Like c-Si, noncrystalline amorphous silicon (a-Si) is semiconducting at 
low pressure and maintains a coordination number A^^ ~ 4. At higher pressures, amor- 
phous silicon becomes metallic and more highly coordinated. Thus, both c- and a-Si exhibit 
pressure- induced polymorphism. 

Liquid silicon (1-Si) is metallic at high temperature, although it possesses a small popula- 
tion of covalent bonds P, 5]. The coordination number ~ 7.3 (for liquid and amorphous 
structures is the average number of neighbors within the first peak of the radial distri- 
bution function g{r)) is lower than occurs in ordinary metals such as aluminum or copper, 
but higher than in c-Si. As expected from its higher coordination number, 1-Si has higher 
density than c-Si, a property that silicon shares with water 
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The structure of supercooled 1-Si and its transition from a metallic dense-packed struc- 
ture at high temperature to a semiconducting open network at low temperature remains 
uncertain. The dissimilarity of the amorphous and liquid states raises the possibility of 
intermediate phases p]. By analogy with the liquid-liquid phase transition (LLPT) in wa- 
ter j^, Q|, where tetrahedral order is also present in the crystalline solid phase, researchers 
speculate [3] such a transition might occur in silicon, between a tetracoordinated low density 
liquid (LDL) and a more highly coordinated high density liquid (HDL). Since the LLPT pre- 
sumably occurs in the supercooled liquid, both LDL and HDL are metastable states, if they 
exist at all. The supercooled HDL state is the metastable extension of the high temperature 
equilibrium liquid (HTEL) state, while the glassy a-Si state is the extension of LDL to low 
temperature. 

Experimental evidence for an LLPT is not conclusive. Pressure induced transformation 
of a-Si [Sj, from a low-density amorphous (LDA) state at low pressure to a high-density 
amorphous state (HDA), similar to that of amorphous ice Q|, suggests the possibility of an 
LLPT in silicon similar to that in supercooled water. Plastic deformation 10| studies indicate 
a transition from a-Si to a fluid state, tentatively identified as LDL, at a temperature around 
lOOOK (the mdfug po„. . T„. ^ 1687/0- Co„ta,„erless supercooling expen,„e,.s 
121 . 1131 . 1141 . |15I | lead to contradictory results regarding the temperature dependence of the 
supercooled liquid properties. 

Computer simulations reveal complete structural detail not available by ordinary exper- 
imental means. However, they suffer from uncertainty arising from approximate models of 
interatomic forces and energetics, as well as limitations of sample size and simulation time. 
To capture the nature of tetrahedral bonding, the Stillinger- Weber (S-W) potential pJI] 
includes angle dependent interactions. This potential shows clear evidence of a first-order 
HDL to LDL transition, in the deeply supercooled regime (with a transition temperature 
Tc ~ 1060/^) Electronic structure calculations on these simulated configurations re- 



port a metal to semimetal transition across the LLPT transition temperature pj|. Other 
empirical_ 
tential 



potentials, such as the Keating potential [19| and the environment dependent po- 
20], find contradictory results, so the nature and existence of an LLPT depends on 
the form of the interaction model [2l| . 

Since neither experiments nor simulations with empirical potentials resolve the existence 
or character of an LLPT in Si with certainty, further investigation using first-principles 
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methods is warranted. The pioneering Car-Parrinello simulations Q] studied sohd c- 
Si and a-Si in addition to high temperature 1-Si. A later simulation of the supercooled 
liquid 23|] in a small = 64 atom cell found increased tetrahedral order and a sudden drop 
in Nc at T=1100K. Our study carries out a more extensive investigation and finds van der 
Waals loops in the pressure- volume isotherms (see Fig. [1]), a strong indication of an LLPT, 
with a critical point around ~1232K. 

Our simulations use the Vienna Ab-initio Simulation Package (VASP) 2S, |25|] together 



with Projector Augmented Wave (PAW) potentials [26|, |27|] in the Generalized Gradient 
Approximation (GGA). Electronic structure calculations of energies, pressures and forces 
used the P fc-point only and were conducted at the default energy cutoff of 245eV. Molecular 
dynamics utilized a 2fs time step, and velocities were rescaled every time step to maintain 
constant temperature. Our molecular dynamics calculations extended to 3ps or more and 
required 2 weeks of computer time for each temperature and volume. 

To test for finite size effects we evaluated the radial distribution function g{r) at system 
sizes A^=100, 200 and 300 and found adequate convergence at = 200 atoms. We trans- 
formed g{r) to obtain the structure factor S{q), using the Baxter method [28| (see Fig. El). 
The peak and shoulder positions agree well with experiment {29!, but the experimental S{q) 
exhibits a slight excess at low q arising from a background effect, causing the experimen- 
tal data to violate a sum rule associated with short distance correlations while the 
simulation data obey the sum rule. 

Ideally a phase transition should be observed as a function of temperature T and pressure 
P. However, the VASP program is restricted to fixed volume simulation, so we monitored 
the reported "external pressure" (derivative of potential energy with respect to volume [30I ) 
and added the kinetic energy contribution pksT to obtain the thermodynamic pressure. We 
also add a constant offset of Ppuiay=^-^ kilobar to account for the Pulay stress, which we 
found to be reasonably independent of temperature and volume. 

At least 9 temperatures were studied at each volume, uniformly spaced between a low of 
T=982K and a high of 1382K, with additional intermediate temperatures included in cases 
where the mean internal energy varied rapidly with temperature. Each starting configuration 
for a given V and T was generated by a lengthy MD simulation prior to the beginning 
of data collection. To improve the efficiency of configuration sampling, we employed a 



replica exchange method 28|, l3l|, l32| to swap configurations between temperatures. Briefly, 
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pairs of configurations at common volume but differing in energy by AE and differing in 
inverse temperature (/3 = l/ksT) by A/3, were swapped with probability exp {ApAE). We 
attempted swaps every lOOfs. 

Having energy and a series of nearby temperatures allows us to 



employ the "multiple histogram" technique 32|, ISSj. Histograms of internal energy at 



fixed temperature Ht{E) are converted into configurational densities of states Q{E) = 
Ht{E)) / C^j, e^^'^'^^~^^^''^'^) where the sums run over the discrete temperatures T at 
which simulations were performed. Consistency between temperatures is enforced by their 
free energies E{T) = -kBT\nZ{T) where the partition function Z{T) = f fi(E)e-^/'=s^dE. 
We thus reconstruct the free energy and its derivative quantities such as pressure and heat 
capacity as analytic functions of temperature. For example the pressure is obtained from 

P(T) = ^ / P{EME)e~^/''-^dE (1) 

where P{E) is the mean pressure observed at energy E. Accuracy of this method depends on 
a thorough sampling of configuration space and on sufficient overlap of the energy histograms 
at adjacent temperatures. 

Results for the pressure are illustrated in Fig. [TJ Error bars are estimates of standard 



error, defined as the RMS fiuctuation in pressure divided by Nind, where we set Nind=W 
as the approximate number of independent samples. Two main features of this plot are: van 
der Waals loops, and negative thermal expansion. 

Van der Waals loops occur when a region of positive slope interrupts the generally negative 
slope of the pressure-volume isotherm. Positive slope of P{V) corresponds to negative 
isothermal compressibility {Kt = — ■^(|^)t)- This thermodynamically unstable state would 
phase separate into high volume (low density) and low volume (high density) phases if 
the system size were sufficiently large. However, negative K-r is permissible in systems of 
finite size owing to the free energy cost of the interface needed to separate the two phases. 
Maxwell's equal area construction can be used to identify the coexisting states below the 
critical temperature, which we identify as Tc ~ 1232K. At T = 1182K, we find coexisting 
states: HDL at density p = 0.053 atoms/A'^ (volume v = 18.9 A^/atom) with coordination 
number A^^ = 5.6; LDL at density p = 0.049 {v = 20.4) with coordination number = 4.2. 
According to the isotherms, HDL is the extension of the equilibrium liquid. 

Negative thermal expansion {ap = ■^(|^)p) occurs when the pressure is a decreasing 
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function of temperature and can be seen at volumes above 18 A^/atom. Negative thermal 
expansion is thermodynamically permissible, and is seen in water close to freezing, for exam- 
ple, but is quite rare. Thermal expansion is related to entropy via (ff-)T = Since LDL 
exhibits negative ap, we see that entropy decreases (counterintuitively) as volume increases. 
Presumably this reflects the higher degree of orientational order in the open tetracoordinated 



network as compared with the disordered high density liquid state 34| . 

Fig. [2] (inset) shows the distribution of the tetrahedral order parameter [sl, 

?* = 1-^ E(cos% + ^)^ (2) 

i<j=i 

where 9ij is the angle formed by an atom with its i*^ and j*^ nearest neighbors. A value 
close to 1 indicates tetrahedral bonding. We evaluate this distribution in the HDL and LDL 
states at temperature T=1182K. At lower densities, the favored LDL structure is indeed 
much more tetrahedral than the higher density HDL or HTEL structures. 

Because the LLPT occurs between metastable liquid phases, below the equilibrium freez- 
ing transition, a crucial question is if our coexisting phases at low temperature are both 
truly liquid. Indeed, since LDL exhibits fairly strong peaks in g{r) and S{q) as well as a 
high degree of tetrahedral order, a central concern is if this LDL phase has not actually crys- 
tallized. Likely candidates for the crystal structure are the usual c-Si diamond structure of 
Pearson type cFS, or the Pearson type cI16 crystal that arises upon recovery at atmospheric 
pressure following high pressure treatment. 

To make direct comparisons we carried out solid-state molecular dynamics simulations. 
For cF8 we took A^=216 atoms at f=20.4 A^, and for cI16 we took A^=192 atoms at 
f=18.5 A'^, the densities being chosen taking thermal expansion into account. Compar- 
isons of g{r) and S{q) are shown in Fig. [31 In order to smear out the Bragg peaks in S{q) 
we smoothly truncated the crystalline g{r) prior to Fourier transformation, so the figures 
understate the strength of peaks in S{q) for the crystalline structures. At T=1182K both 
crystal structures remain far better ordered than the LDL. Additionally, for both g{r) and 
S{q) there is no systematic resemblance between LDL and either crystal structure. 

Comparing qt distributions at T=1182K (not shown) we found the distribution for cF8 
was more sharply peaked close to qt=^ than that for LDL, but the distribution for cI16 
was rather similar to that of LDL. Coordination numbers also were similar between LDL 
and cI16. We examined the statistics of atom rings. For the ideal crystal structures, both 
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cF8 and cI16 have only even length rings, and no rings of length less than 6. For purposes 
of comparison we take the ratio of the number of rings with length 5 or 7 to the number 
with length 6 or 8. Values of (iVg + N7)/{N(i + Ng) at T=1182K were: 0.00 for cF8, 0.08 
for cI16 and 0.50 for LDL. We also found some rings of length 3 and 4 in LDL, but not in 
either cF8 or cI16. Finally, we visually inspected atomic configurations and saw no hints of 
crystallinity in LDL 

The electronic density of states (DOS) governs important material properties including 
electric conductivity. Fig. H] illustrates the DOS of the equilibrium HTEL, metastable LDL 
and HDL, and c-Si states. Liquid state DOS are presented as the average of five independent 
configurations. The Fermi energies Ep of HTEL and HDL are similar because their atomic 
volumes are similar. Likewise the Ep^s are similar for LDL and cF8, but both lie below those 
of HTEL and HDL because their atomic volumes are greater than for HDL and HTEL. The 
bottom of the LDL valence band (around -7eV) lies slightly below cF8, reflecting the slightly 
higher coordination of LDL. HTEL and HDL have high densities of states at Ep impl ying 
metalhc character. LDL has a deep pseudogap at Ep suggesting semimetallic character |l8 |. 
c-Si exhibits the characteristic bandgap of a semiconductor. 

In conclusion, we find that deeply supercooled liquid silicon undergoes a liquid-liquid 
phase transition, separating into a high density highly coordinated metallic liquid and a 
low density low coordinated semimetallic liquid. Our calculation reveals structural detail 
not currently available from experiment in this temperature range. Because we employ 
first -principles forces, we have confidence in the validity of the model. The critical point 



lies close to T = 1232K, not far from estimates from experiment [8|, llO| and simulations 
based on empirical potentials . HTEL and HDL share characteristics such as low atomic 
volume, high coordination, moderate qt and metallic conduction, consistent with the notion 
that HDL is the metastable extension of HTEL into the supercooled regime. Likewise c-Si 
and LDL share characteristics such as high atomic volume, low coordination, large qt and 
low DOS at Ep, consistent with LDL exhibiting local tetrahedral order similar to c-Si. 
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FIG. 1: (color online) Pressure- volume isotherms of liquid Si. Data points are calculated using 
Eq.E 
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FIG. 2: (color online) Liquid structure factors S{q), radial distribution functions g{r) and tetra- 
hedral order parameter {qt) distributions. High temperature simulation at temperature T=1815K 
and volume per atom v = 18.2 A'^ is compared with experiment 29]. Low temperature simulation 
at T=1182K shows coexisting structures, HDL at v=l8.9 and LDL at v=20A A^. 
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FIG. 3: (color online) Comparison of low-density liquid structure with competing crystal structures 
at T = 1182K. 




FIG. 4: (color online) Electronic densities of states (states/eV/atom). Vertical dashed lines locate 
Ep. 



